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A SPECTRAL METHOD FOR HALF-INTEGER SPIN FIELDS 
BASED ON SPIN-WEIGHTED SPHERICAL HARMONICS 

FLORIAN BEYER, BORIS DASZUTA, AND JORG FRAUENDIENER 


Abstract. We present a new spectral scheme for analysing functions of half¬ 
integer spin-weight on the 2-sphere and demonstrate the stability and conver¬ 
gence properties of our implementation. The dynamical evolution of the Dirac 
equation on a manifold with spatial topology of 8 2 with a pseudo-spectral 
method is also demonstrated. 


1. Introduction 

In a recent paper [3] we presented a spectral method for tensorial partial dif¬ 
ferential equations on geometries with a spherical component. We showed how to 
implement the ‘eth’-formalism based on spin-weighted spherical harmonics follow¬ 
ing the work by Newman, Penrose and others [9,23]. The restriction to tensorial 
equations implied the use of spin-weighted functions with integer spin-weights. In 
the present work we extend our method to also include functions with half-integer 
spin-weight on spherical geometries. In particular, we show how the evolution of 
dynamical half-integer spin fields can be accomplished. Our motivation for doing 
so stems mainly from General Relativity: spinorial fields describing elementary 
particles such as electrons or neutrinos are fundamental sources for the Einstein 
equations and are studied for instance as test fields on Kerr-Newman backgrounds, 
see [10] and references therein. However, there are also interesting developments 
in condensed matter physics in relation to the description of graphene [5, 39] for 
which our method might be relevant. The numerical treatment of the Dirac field is 
particularly interesting as its evolution can be rather peculiar and sometimes even 
counter-intuitive as visualised by Thaller [34,35]. 

It is well known that when working with § 2 (and other compact geometries) 
coordinate singularities can lead to instabilities that spoil the accuracy of numerical 
schemes. These issues can be avoided by working with spectral methods where 
coordinate-independent manipulation of expressions in terms of a well-defined basis 
of functions upon which the action of differential operators reduces to algebraic 
manipulation is possible [4,11]. Thus, issues of coordinate singularities are dealt 
with automatically by the method. 

Our spectral method makes use of the spin-weighted spherical harmonics (SWSH) 
[9, 23, 25, 26] which can be shown to be equivalent to the well-known Wigner D- 
functions [6,29,33] and thus provide a complete, orthonormal basis for L 2 (SU{ 2)). 
In the description of half-integer spin fields, we have found it convenient to work 
directly with the half-integer SWSH due to the particularly simple action of the 
associated <5 and <5' differential operators which may be thought of as constituting 
covariant derivative operators on the sphere [9]. As the action of the aforemen¬ 
tioned operators on the SWSH reduces to simple algebraic manipulation we may 
translate PDE systems into coupled (infinite dimensional) ODE systems and thus 
construct a spectral evolution scheme for an initial value problem (IVP) of interest. 
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Our work extends that of [3,15] into a spectral algorithm for spin-weighted spher¬ 
ical harmonics (SWSH) with half-integer spin-weight s. The method we present in¬ 
herits several desirable properties from [3,15]. In particular it is theoretically exact 
if a minimum number of grid points are used at a given band-limit L. Furthermore 
the same algorithmic complexity of 0(L 3 ) is achievable for spectral transforma¬ 
tions. In brief, [3,15] seek to calculate values for the integer SWSH over § 2 by 
mapping the sphere into the 2-torus (thus allowing for Fast Fourier Transforms to 
be used) and relate the SWSH to the reduced Wigner d-functions evaluated at tt/2. 
For a short summary of past work related to the integer case see [3]. The approach 
to the half-integer case is in principle similar, however suitable modification of the 
2-torus map must be made in order to account for the periodicity of spinor fields. 
In addition, the recursion relation [36] that allows for evaluation of the Wigner 
d-functions must be modified. 

As a toy model we numerically explore the dynamics associated with a Dirac 
equation on a 2-dimensional manifold with spatial topology of § 2 as an IVP. This 
will allow for a test of the spin-weighted spectral transformations we present in 
the context of evolution equations. Common techniques for treating the numerical 
problem of the Dirac equation consist of: FD schemes formulated on a flat-lattice 
in configuration space [12, 13], on a grid within a finite-volume in momentum- 
space [22] and using methods based on the split-step operator technique [7, 8, 20, 
21], Particular to the FD approach special care must be taken so as to avoid the 
Fermion-doubling problem [24]. Elimination of spurious modes introduced to the 
solution may be accomplished by means of nonlocal approximation for the spatial 
derivative operator [32,38] or by staggered-grid schemes [12,13]. As remarked upon 
in [13] the issue of spurious modes is not particular to the Dirac equation but 
can occur whenever a symmetric FD approximant is used for a first derivative on 
a uniform grid. We seek to avoid the above issue entirely by making use of the 
global approximation to functions and derivative operators that spectral methods 
provide [14]. 

The equations of motion (EOM) that we derive for the Dirac equation in the 
aforementioned geometric setting result in a coupled system of variable-coefficient 
linear hyperbolic PDEs. When performing a SWSH decomposition this results 
in product terms that must be further simplified. From past experience [3] we 
have found that while it is possible to make use of Clebsch-Gordan expansions (see 
sec. 2), i.e., working entirely in the space of coefficients, it is far more convenient 
to instead work with a pseudo-spectral method (see sec. 5.1), and thus we follow 
the latter approach in this work. 

This paper is structured as follows: In sec. 2 we recall basic properties of the 
half-integer SWSH together with the (3-formalism. In sec. 2.1 and 2.2 respectively 
we discuss the forward and backward spectral transformations of spin-weighted 
functions on the sphere. In sec. 3 we present the required modification to the re¬ 
cursion relation for computing the reduced Wigner functions at tt/2. Subsequently, 
in sec. 4 we perform consistency checks on our implementation of the algorithm. 
Error pairs and the property of exponential convergence are analysed. In sec. 5 we 
proceed by numerically solving the 2 + 1 dimensional Dirac equation as an IVP on 
a curved geometry, specifying to cases with spatial topology of § 2 . We find EOM 
adapted to the spin-weighted formalism. In sec. 5.1 we briefly recall the pseudo- 
spectral method. In sec. 5.2 we construct numerical solutions to the EOM and 
inspect convergence properties of the numerical solutions obtained together with 
conserved currents. Sec. 6 concludes. 
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2. SWSH SUMMARY 


In this section we briefly summarise key properties of the spin-weighted spherical 
harmonics (SWSH) that we will make use of in the presentation of our half-integer 
spectral algorithm and in the pseudo-spectral method for solution of dynamical 
systems. For further details we refer the reader to [3,6,9,23,25,26]. 

Geometrically, the spin-weighted functions are sections of certain line bundles 
over the 2-splrere 8 2 . They have a representation in terms of ordinary functions 
over patches of § 2 which depend on the choice of coordinates and the choice of an 
orthonormal frame at points of the patch. The behaviour of this representation 
under changes of the frame is captured by the spin-weight assigned to the function. 
Here, we will deal only with representations of the spin-weighted functions given in 
terms of the standard polar coordinates (i9, <p) and the orthonormal frame defined 
in terms of the coordinate derivative vectors 

— — -cU. (2.1) 

sin tit 

Equivalently, these vectors can be obtained as real and imaginary parts of the 
complex linear combination 


m = 


(e$ ie v ) 



i 

sintf 



which, together with its complex conjugate, satisfies the orthonormality relations 


m • m = 0, m ■ m = 0, m • m = 1 

at all points of § 2 covered by the polar coordinates with respect to the standard 
metric on § 2 . 

Every smooth function s f on § 2 with spin-weight s can be expressed as: 

l i 

sf(ti),ip)= lim E V sdlmsYlmi $,¥>), (2.2) 

L—too z — J z " 
l—\s\ m——l 

where the s Yi m are the SWSH, which form a complete orthonormal basis for the 
complex vector space of spin-weight s functions. 

Recall that the action of the 3 and 3' differential operators on s f serves to raise 
and lower the spin-weight respectively; when viewed as maps with respect to polar 
coordinates we have: 

3 : s / —t s +if, 3 [ s /] = d$[ s f] - icscdc^s/] — scotti) s f, (2.3) 

3' : s f -)• s _i/, 3' [ s /] = dfilsf] + icscti)dip[ s f] + scotd s /. (2.4) 

The SWSH may be written explicitly as: 

s Y) m (tf, ip) = y^l'expCim^) d l sm (i9), (2.5) 

where d( rin (i?) is the reduced Wigner d-function: 


min(Z+m,Z— n) 

d l mn w = E (-iy- m+n 

r—max(0,m—n) 


yj(l + m)\(l — m)\(l + n)\(l — n)\' 

r\(l + m — r)\(l — r ~ n)\(r — m. + n)\ 

i9 


X cos 


2l—2r+m—n 


2 1 sm ' 


2 r-m-\-n 


Under complex conjugation the SWSH satisfy: 


( 2 . 6 ) 


Y lm = _ s Yi,_ m 




(2.7) 
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The SWSH satisfy the orthonormality relation: 

sYl^mf) / / 3 ^imi (^) l 2 m2 V 3 ) ^ dif di 1 i 2 S mirrl2 . 

Jo Jo 

( 2 . 8 ) 

A particularly useful property the SWSH possess is that under the action of <3 and 
3' the action of Eq. (2.3) and Eq. (2.4) reduces to a ’’ladder algebra”: 

3 [ s Yi m ($,</>)] = -yj{l- s)(l + s + l) s+:L Y) m (tf, ip), (2.9) 

3' [ s Yim{fi, V 5 )] = \/(l + s)(l - s + 1) ' s _ibm(l?, <f), (2-10) 

which we will exploit when performing decompositions of dynamical equations, 
allowing for a map from a PDE system to an (infinite dimensional) ODE system 
(see sec. 5). In addition to this we have the commutator expression: 

[3,3'] a Y lm (0,(p) = -2a a Y lm (0,tp). (2.11) 

A product of two SWSH with spin-weights Si and S 2 is a function of spin-weight 
si + S 2 and, therefore, it can be written as a finite linear combination of SWSH 

(A < / 3 )s2-^j 2,m 2 (d, v ) = Ci(si,li,mi; s 2 ,Z 2 ,m 2 )( Sl+S2 )l),( TOl+TO2 )('i?, <p), 
le A' 

( 2 . 12 ) 

where A' := {max(|Zi — Z 2 1, |si + S 2 1, \fn\ + TO 2 I), ..., h + I 2 } and the coefficients 
are related to the usual Clebsch-Gordan coefficients, see e.g. [3,29]. 


2.1. Forward transformation. We now describe our numerical algorithm for 
evaluation of the forward transform T : s f 1 — > (sUim)- As a first step we intro¬ 
duce the notation A l mn := d l mn (n/2) which allows for the rewriting of Eq. (2.6) 
as [28]: 

1 

dLnW = i m ~ n J2 A l qrn e~^A l qn , (2.13) 

q——l 

following from a factoring of rotations [36]. In particular, note that d l mn {’d) = 
dlnnid) (cf. Eq. (2.6)). We defer the details of how the A elements are calculated 
together with their symmetry properties to sec. 3. 

Define the functional: 


/• 2"7T /» 7T 

Imn[sf(d,ip)]:= / / e~ lm ^e~ m,p s f(d, ip) sin ddtfdip. (2.14) 

Jo Jo 

This integral can be evaluated exactly - we now describe our method which is based 
on [15]. Combining eqs. (2.2), (2.5), (2.13) and (2.14) results in: 


s dim — 1 


2Z + 1 
47T 


V A 1 I A 1 


21 + 1 
47T 


q=—i 
1 1 


A qs J qrn A 


1 

qm 1 


9 = 1/2 


(2.15) 

(2.16) 


where 

jqm := I qm + (-1 ) 2l + S+m I-q,m. (2.17) 

In order to make use of existing FFT (fast Fourier transform) algorithms in the 
evaluation of the quadrature in Eq. (2.14), we must extend a spin-weighted function 
s f with data sampled on the domain D := {($, (/>)|$ £ [0,7t],</j £ [0, 27t]} to a 47 t 
periodic function on D := {($,</>)[$ £ [0,47 x),ip £ [0,47t)}. This is accomplished 
using the intrinsic symmetries of the sYimid, ip) together with the expansion ansatz 
of Eq. (2.2). For convenience we define the domains D\, .. ., Dym (see Fig. 1). 
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We have D = D\ upon which s f{ft, ip) = s fi{ft,<p). For Du put s /n(i?, ip) = 


Di 

D„ 


,fn 

D m 

Div 

sflll 

,fiv 

D v 

Dyi 

,fv 

si VI 

Dvn 

Dym 

s^VII 

sfvm 



Figure 1. Domains for function extension. The shaded region 
corresponds to data that are to be extended to unshaded regions. 


— s fi{ft, p — 2 tt). Next, define s g : D\ U Du —> C by 

^’ v)€Dl . 

v ‘ \, In («,*>) («,¥>) €i>ii 

Now, define s h : D\ U D\\ U Dui U Div -A C by 

^ , = (i?, ip) £ Di U Du 

\ (~l) s+1 s g{2n -ft,(p- 7r)mod47r) {ft, ip) £ Du\ U Diy 

Finally, we arrive at the extended function S F : D —» C defined by 

p,^ ,_f s h{ft, <p) (ft, ip) £ D\ U Du U Dm U Di V 

^ | —shift — 27r, <^) ($, </?) £ Dy U Dyi U Dyn U Dym 

The smooth function s F{ft, p) is 47 t periodic in ft and ip, thus we may write: 


(2.18) 


(2.19) 


( 2 . 20 ) 


K M 


s F(ft,ip) = EE s-^km^ 




( 2 . 21 ) 


k —0 m =0 


In order to evaluate (2.14) we now identify 

/*27T /»7T 

Imn [af{ft, p)\ = Imn [ s F(ft, p)] = / e~ irn ^ /2 e~ mv/2 s F(ft, tp) sin I? dftdip, 

J o do 

( 2 . 22 ) 

where we can extend the indices to all integers, i.e., m,n £ Z. This is permissible 
as both functions are equal in the region of integration. Upon substitution with 
Eq. (2.21) we find: 


with 


Imn [ s F(ft, ip)] = ^2 sFij W# (i - m)W ip (j - n), 
i,j=0 


W#(t) 


F T ®! 2 sin ftdft 


±i- 
ql+T 
4—T‘ 


T = ± 2 
r £ Z \ {±2} ’ 



27T 


P = 0 

P£ Z\{0} ' 


(2.23) 

(2.24) 

(2.25) 
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Equation (2.23) amounts to a two-fold discrete convolution in spectral space. By 
the convolution theorem, this implies that we may equivalently consider pointwise 
multiplication of the inverse transforms of W# and W v with S F. Let N# and N v 
denote the number of grid points over $ and ip on the containerization of S 2 that 
the function s /(i9, ip) is to be sampled at. Upon extension to the doubly 47T periodic 
domain (see Fig. 1) we work instead with the extended function S F{ $, ip). Further¬ 
more, in order to simplify numerical construction of the extension, the number of 
samples on the extended domain is taken to be N# = 4 N# — 3 and N' = 2N ip — 1. 
Correspondingly, the spatial sampling intervals are now given by At? = t and 
Atp = or der to satisfy the Nyquist condition, such that spurious aliasing 

does not occur, we impose N# = N v = 2 (L + 2— 1/2) 4-1, where L is the harmonic 
that the function s f{&,p) is band-limited to 1 . 

Under these choices a convenient form of the quantities discussed above that is 
directly amenable to numerical work is given by: 


W*(7) ■•= 

W v (s) := 
together with: 


v i9 — z 

E 

5=0 



N f 2 / N' — 1 

E e^ l2 wJ^~ -C 

C=o V 


7 =0,... N# — 2 

(2.26) 

£=o,...iv;-2 

(2.27) 


s F(t?,<p)] = 


1 


N'-2N'-2 


E E 


Eo E=0 

x W#('y)W ip (e) a F('y,e). 


e -ifi'yA r d/2 e ~ivsAip/2 


(2.28) 


where p £ 0,..., N# — 2, v £ {0,..., N' — 2}. Note that in this approach, we 
must discard all values of I^ v for which p and v are even integers. In order to 
recover I mn from Eq.(2.28) we take p = |( N# — Am — 1) and v = \{N' —An — 1). 
Overall we find an algorithmic complexity of 0(L 3 ) as the integrals may be 
evaluated exactly in 0(L 2 logL) operations by performing a two-dimensional FFT 
and each component of A/ m (as required by Eq. (2.16)) can be computed using 
0(1) floating point evaluations - see sec. 3 and also [15,36]. We remark that if the 
analysis of strictly real data is desired then it is possible to attain a linear increase 
in the execution speed of transformations (specifically the FFT component by a 
factor of approximately 2). However, as we are primarily interested in applying 
transformations for the solution of equations of motions that govern the dynamics 
of complex fields we do not explore this further. 


2.2. Backward transformation. We now describe the algorithm for evaluation 
of the backward (inverse) transform J 7-1 : ( s a; m ) i-A s f. The backward spherical 
harmonic transform maps the expansion coefficients s a; m , for |s| < l < L, to a 
function on (a dense subset of) § 2 . Because we are working with band-limited 
functions we can, at the analytical level, perfectly reconstruct the original function. 
To this end, Eq. (2.2) must be evaluated. As the inverse transform does not contain 


1 More efficient samplings on the sphere may be possible [19] however the overall asymptotic 

complexity of our proposed algorithm will not be affected and thus we do not explore this issue 

here. 
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integrals, issues of quadrature accuracy do not arise. Define: 

L 


Kmn [sain] ~ R _ " 


l=\s\ 


477 


which allows for: 

,M<P)= E E e imS e inip K mn . 

Use of Eq. (3.1) permits a rewriting of Eq. (2.30) as: 


Kmn [s^Zn] — 


i S " £j=|«l V A -m,ssain Al m , n ™, > \ 

i s ” n Ef =k , JW^ l ) 2 l+S+n ^ l m,ssain^ m , n rn<-l 


(2.29) 


(2.30) 


, (2.31) 


which allows for a reduction in computation time. In Eq. (2.30) we may use an 
FFT directly, discarding values of s /($, tp) for which $ > 77 and <p > 2n. If the 
input data to the FFT library is Hermitean then another linear increase (~ 2) in 
the execution speed of a transform is possible however, this again does not change 
the overall algorithmic complexity 0(L 3 ). 


3. WlGNER A AND RECURSION 


Here we answer the question of how to compute the A l mn required in the forward 
and backward transforms. We will base our numerical algorithm for computation 
of an arbitrary A l mn on recursion and exploitation of symmetries. 

First we note the symmetry properties (inherited from dr mn ($)): 


(—l) l+n A l mn 

(3.1) 

(-l) l ~ m A l mn 

(3.2) 

(-ir~ m A i nm . 

(3.3) 


Using an approach similar to [36] one can derive an analogous Trapani-Navaza (TN) 
style recursion from Eq. (2.6) for half-integer l,m,n values: 


A l u = 2 - 


A' = 


1 ( 21 - 1 ) 


A 


1-1 


2(1 + m)(l + to — 1) 1 




A' = 

mn 


2 TO 


\/(l — n)(l + n + 1) 


A 


1 

m,n+1 


l (l-n-l)(l + n + 2) , 

(l-n)(l + n+ 1 ) m ’ n+2 


A‘ 

I — 


21 


A‘ _i 

2 / + 1 


(n = l — 1) 


< 1 - 2 ) 


(3.4) 

(3.5) 

(3.6) 

(3.7) 


(3.8) 


Note that by making use of the symmetries provided by eqs. (3.1)-(3.3) similar re¬ 
cursion relations may be constructed connecting different combinations of subscript 
indices. 

We visualise the possible values of A l mn up to some maximal band-limit L as 
being arranged in a square pyramidal lattice with A l mn values corresponding to 
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l = 1/2 and —1/2 < m,n < 1/2 occupying the top-most plane; l = 3/2 and 

—3/2 < m,n < 3/2 the next plane down and so forth. In a recursive approach one 

1 

can thus initialise with a single value of A mn and with eqs. (3.4), (3.5) and (3.8) 
(together with symmetries) compute those values of A l mn constrained to the surface 
of the pyramidal lattice. For each fixed l value those A l mn that occupy the interior 
of the pyramidal structure can be computed using Eq.(3.6) and Eq.(3.7) (together 
with symmetries). 

We note that symmetries allow for a reduction in the total number of elements 
that must be calculated explicitly via recursion to ((2/ + l)/2) 2 entries for a fixed l 
plane and y^( 2L + 3)(2L + 1)(L + 1) total entries for a choice of maximal band-limit 
L. 

In our implementation we initialise with -x/i ~ l/v^ and iterate such that 
for a given l the to, n indices obey the condition (n < 1/2) A (|m| > |n|) A (to > 
-!/2). 

Working with double precision arithmetic, we have found that the above scheme 
remains stable up to a band-limit of L ss 5173/2. Due to exponential convergence [4] 
this L will in practical situations be far above the resolution required to accurately 
sample fields for evolution equations and thus is not a concern for this work. We 
note however, that by instead working with ratios such as r A(, n = A l mn /A l m _ ln 
we have verified that it is possible to construct a stable type of hybrid recursion in 
analogy to [3], for L > 5173/2. 


4. SWSH CONSISTENCY CHECKS 

As a first check on the consistency of the half-integer SWSH algorithm presented 
above we construct error pairs. Here one populates coefficients with random data 
whereupon a transformation is applied so as to construct the spatial representa¬ 
tion of the corresponding function; this spatial function is transformed back to 
coefficients and the associated error may be examined. This procedure may be 
summarised as: 


s^lr 


s /($,¥>) 


s^lr 


(4.1) 


The real and imaginary parts of s a/ m we generate by sampling from the continuous 
uniform random distribution on the interval [—1,1). The numerical error associated 
with the specification of Eq. (4.1) is shown in Fig. 2(a) and Fig. 2(b). We see that 
our implementation is indeed consistent, accurate and stable. Furthermore the 
scaling L(n) of e re irms (Fig. 2(b)) for the half-integer spin-weight uniform random 
data matches the scaling observed for integer spin-weight Gaussian random data 
observed in [15]. 

The product of two spin-weighted functions Sl f and S2 g on § 2 is a function with 
spin-weight si + S 2 - If Sl f and S2 g are comprised of a finite number of SWSH then 
their product is also and thus by selecting a sufficiently large L we exactly sample 
the resulting product function. We can also check the property of exponential 
convergence; to this end define Ai := ( s ai m ) m = \s a im |/(2Z + 1), a measure 
of the average magnitude of coefficients at a fixed l value. We expect that given 
smooth test functions Ai should behave as Ai ~ aexp(— nl) (a, k £ R) for large 
l [4,16]. Smooth half-integer spin-weighted functions may be constructed by taking 
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(a) (b) 

Figure 2. (Colour online) Averages of numerical error pairs 
when performing the transform procedure of Eq. (4.1). A fixed 
spin-weight s is chosen and 5 sets of s aim are randomly generated as 
described in the text. In both sub-figures the band-limit of trans¬ 
forms is given by L(n) = (2 n — 1)/2. (a) Averaged maximum, abso¬ 
lute relative error e re i := (rnax; m |1 — s^lm/s&lml) ■ The thick black 
line corresponds to a scaling of L(n) 2 . (b) Averaged maximum, ab¬ 
solute relative rms error e re irms := (( s ai m - s d; m ) rms /( s d; m ) rms ). 
The thick black line corresponds to a scaling of L(n). In (a) and 
(b): Green ”+” corresponds to s = 3/2; blue ”o” corresponds to 
s = 13/2; black corresponds to s = —1/2. 


a finite number of SWSH and modulating by the exponential of a smooth spin- 
weight 0 function. Introduce: 


o = i 0^,-1 (*?,¥>) + IT 0 ^ 1 ( 1 ?,^) (4.2) 

0 gift, f) = exp (- 0 g(fl, <f)) (4.3) 

1/2 f{^,f) = 1-3 1 / 213 / 2 , 1 / 2 ( 1 ?, f) + i 1 / 2 1 3 / 2 ,- 1 / 2 ( 1 ?, f) (4.4) 

1/2 /($, f) = 0 g($, f) 1/2 f(&, f) (4.5) 

3/2 h{&,(f)= [ 1 / 2 ^ 11 / 2 , 3 / 2 ( 1 ?, V^)] 3 (4-6) 

_i/ 2 fc('i?, f) = 0.7 i —i/2^5/2,— 1 / 2 ( 1 ?, f) + 0.9 - 1 / 2 I 3 / 2 , 1/2 ($, f) (4-7) 
- 1/2 fc(i?, f) = o<?(^ f) - 1/2 fc(i?, f) (4.8) 


In Fig. 3 we show Ai for various combinations that test the various properties 
outlined. In particular, we observe that functions comprised of a finite number of 
SWSH are completely captured if L is greater than the highest harmonic of the 
linear combination. Products of spin-weighted functions (each individual function 
comprised of a finite linear combination of SWSH) also are exactly resolved to 
within numerical error for appropriate L. Furthermore, we observe the crucial 
property of exponential convergence for smooth functions. 
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Figure 3. (Colour online) The quantity Ai is depicted and corre¬ 
sponds to: blue ”o”, 1/2 f(fi,<p)\ red 1 / 2 f('&, <p) og{$, <f)\ black 

1/2 /(tf,¥>); green 3/2 cyan _i/ 2 fc(tf,/?). All 
transforms are performed at a band-limit L = 129/2 and s = 

1/2, 1/2, 1/2, 3/2, —1/2 respectively. The functions i/ 2 /($, ¥>), 
i/ 2 /($, V 3 ) V 7 ) and 3/2 h{$,ip) are completely sampled as the 

average magnitude of coefficients ( A [) falls to numerical round-off 
for l well below L. For i/ 2 /(i?, <p) and -\/ 2 k ;($,/?) we observe the 
expected exponential decay in Ai. (See text for details) 

5. 2 + 1 Dirac equation 

In order to test the SWSH half-integer algorithm we now numerically solve the 
Dirac equation on a 3-dimensional Lorentz manifold (•/#,</) with topology lx§ 2 . 
Since the Dirac equation on a 3-dimensional manifold is less familiar than its 4- 
dimensional counterpart we present a detailed derivation in Appendix A and simply 
state the result here. 

We consider the manifold ~KxS 2 with the metric 

g = dt®dt — ■$, <p) (di? ® d$ + sin 2 fidip ® d/s) , (5.1) 

where (d, ip) are standard polar coordinates for the 2-sphere. The function PF{t, $, ip) 
is a conformal factor relating the induced metric at every instant of time t on S 2 
to the standard metric of the unit 2-sphere. 

According to Appendix A the Dirac equation on is given by the following 
two equations for a spinor 1 / = [ip-,ip + ] T 

d t ip- = - ^d'Tip + - •F3V+ + y- 0-, (5.2) 

d t ip + = i pip+ - ip- - Tdip- + yrip+. (5.3) 

The spinor components tp± have spin-weight ±1, while T has vanishing spin-weight. 
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5.1. Pseudo-spectral method. While performing a spectral decomposition of 
variable coefficient PDEs such as Eq. (5.2) and Eq. (5.3) by assuming functions 
may be expanded as in Eq. (2.2) and then products decomposed as in Eq. (2.12) 
is possible (see [3] for examples of such expansions in the integer SWSH case) we 
have found that instead it is far simpler 2 to use a pseudo-spectral (PS) approach. 

We outline this as follows: Given the coefficients Sl /; irrai and s 2 gi 2 m 2 which rep¬ 
resent the functions Sl f and S2 g sampled at a band-limit L then the coefficients 
Sl+S2 ai m corresponding to the associated point-wise product Sl f ■ S2 g can be calcu¬ 
lated by performing the transformations: 

F • ^si/limi j t si/j F '■ (sijlimi) P± s 2 9 i 

subsequently taking the pointwise product and transforming: 

F : Slf s 2 9 P”^ ( s 1 +s 2 0'lm ) J 

we find an approximation to an expansion utilizing Eq. (2.12) directly. We empha¬ 
size that this method also easily allows one to take into account the action of the 
<5, S' operators on spin-weighted functions by embedding their action as multiplica¬ 
tion (see Eq. (2.9) and Eq. (2.10)) in coefficient space, together with taking account 
of their spin raising and lowering properties when transforms are performed. 

We now recast the PDE system of Eq. (5.2) and Eq. (5.3) as an (infinite dimen¬ 
sional) ODE system using the PS method. It is convenient to define the auxiliary 
fields: 


3- :=F&il>+, 

: = ^+S 'F, 4-_ : 

: = *_* 
T 

S +:=F&I>-, 

$+ : = ^ip-dF, ±+ : 

II 

±- 


where ± corresponds to a spin-weight of ±1/2. Using these fields we find that 
Eq. (5.2) and Eq. (5.3) can be written as: 

Ip— ,lm(P) = i/LI/f— Jm {t) — ± ^ — 

„ (5.5) 

' t P+,lm(p) = (f) ^ > +,Im(^) 1=J +,/m(^) ± ^ + 


We note that in the case of F = 1, eqs. (5.2) and (5.3) yield a constant-coefficient 
PDE system which may be directly decomposed as: 


1p-,lm(t) = -i H1p-,lm{t) - ( l + - ) 1p +! lm(t) 


1p+,lm(t) = i g.1p +) lm(t) + [l ± - %p- t lm(t) 


(5.6) 


In this case, a solution is readily arrived at: 
1 


ip±,lm(t) = 


Ul 


(jjJi cos ioJit) ± i/rsin(w;t)) V’±,im(0) ± (I ± —) sin (u>it)ip ,lm (0) 


(5-7) 

where w; := \J(l + ^) 2 ± £t 2 '. 

When F ^ 1 we numerically construct the spatial representation of the function 
using Eq. (2.5) with spatial sampling to coincide with the half-integer SWSH trans¬ 
formation. We will also find it useful to introduce the rescaled current component: 

nt.OM := (5.8) 


2 


Indeed the pseudo-spectral method generalises more readily to nonlinear problems. 
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which with the PS method can be computed using the rescaled fields 'ip± := ip±/E 2 . 
From Eq. (A.11) and Eq. (5.8) the probability Q(£ T ) can be computed via: 

<2( S r) = + i’-Amtp+Jm) ■ ( 5 - 9 ) 

l,m 


5.2. Numerical solutions and convergence tests. In this section we examine 
numerical solutions to the EOM eqs. (5.2) and (5.3) under three conditions: The 
case T = 1 (corresponding to a S 2 spatial geometry), where we numerically solve 
the system of equations given by the decomposition of Eq. (5.6); the case of a 
static deformation (time-independent) of a time-dependent T. In the latter 
two cases we work with the decomposition of Eq. (5.5). In each case we are free 
to select arbitrary initial data for ip±. As convergence tests must be performed 
for any numerical calculation we begin by briefly summarising the procedure we 
follow (based on [2]) for doing so, subsequently presenting our results. Aside from 
presenting convergence for the fields ip± we also verify that the continuity equation 
for the current is obeyed and probability is conserved in each case examined. 

In using the PS approach of sec. 5.1 to compute the time-evolution of ip± we solve 
a truncated (i < L) ODE system for ip±,im{t). Suppose that we use an explicit, 
temporal integrator with a fixed time-step St. The numerically calculated solution 
for this choice of St we denote as St). Assume that there exists a Taylor 

expansion of St) about the exact solution ip±j m {t) with error constant 

terms {£)}: 

OO 

St) = Ip±j m (t) + ^2 St n E n . (5.10) 

n—1 

Suppose now that the integrator is of order p so that {Ej}j <p = 0. Thus: 

ip±,im(t; St) - ip±,im(t) = St p E p + 0(St p+1 ). (5.11) 

By successively rescaling St with a constant (here 2) and comparing ip± ,im(t ; St/2 k ) 
with ip±,im{t) we find: 

2 fcp | ip±,im (t; ^St^j - | —t St p E p . (5.12) 


Equation (5.12) may be used in the case when an analytical solution is known. 
If this is not the case we may instead successively compare St/2 k ) and 

St/2 k+1 ) which yields: 


2 fcp jVht.im (t; ^St^j - ip±,i m } -» (i - St p E p = 8t p E p , 

(5.13) 

this is known as a self-consistent convergence test (SCCT) [2], Note that in both 
cases (Eq. (5.12) and Eq. (5.13)) we may also make comparisons of functions in the 
spatial representation by transforming from coefficient-space. In particular, we will 
consider the maximum, absolute, relative error metric: 


e r (0±(i,i?,y0, V 5 )) : = max 

m,n 


(f>±(t,mA'&,nAip) 
ip±(t,mA’&,nAip) ’ 


(5.14) 


which allows for the comparison of the two solutions <p±, ip± at sampling nodes 
mAi?, nAip. We measure error in this manner such that convergence properties of 
solutions and any potential instability that may result from the SWSH transforma¬ 
tion algorithm may be simultaneously inspected. 

Having described the methods we use to test convergence of solutions we now 
examine some example cases. We fix the mass-parameter to be /i = 1.2 at the outset 
and choose a temporal range of t £ [0, 5] throughout. These choices have been 
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made so as to allow for observation of multiple oscillations of modes (see 

Eq. (5.7)) while simultaneously avoiding exact integer multiple of the frequencies 
Wj/( 2ir). Note that all specified initial conditions ip± | (=0 are constructed with data 
as indicated below and then the corresponding coefficients ip±,im\ t= o rescaled by an 
overall factor so that Q is normalized to 1 at t = 0. For each numerical evolution 
the temporal integrator we choose is the standard, explicit Runge-Kutta 4th order 
method (RK4); thus p = 4 in Eq. (5.12) and Eq. (5.13). 

For J = 1 we select the smooth initial data: 

|f= 0 = ip+it, 1 #, <p)| t=0 = 1/2 f{&, <P) (5-15) 

with _i/ 2 given by Eq. (4.8) and 1 / 2 /($, <p) given by Eq. (4.5). The system 
to solve is thus specified by Eq. (5.6) and Eq. (5.15) upon mapping the latter 
ip±(t, 1 ?, V ? )lt=o to coefficient space using the SWSH transformation algorithm. The 
results of convergence tests (based on Eq. (5.12) with the analytical solution of 
Eq. (5.7)) of our numerical solution are shown in Fig. 4. We find excellent agreement 
with the expected 4th order convergence in time for initial data and parameters 
specified. In particular, we see that the numerical error of the solution is dominated 
by the temporal discretization. In Fig. 5(a) and Fig. 5(b) we provide representative 



k 

Figure 4. (Colour online) Convergence test of numerical solu¬ 
tion ip± (spatial representation) for the case T = 1. £ r ('tp ±) := 
e r St k , 1 ?, ip), ip±(t, •d, </?)) is displayed at t = 5 for a number 

of steps N = 2 fc x 100. i/>±{t,’d,<p) is calculated from Eq. (5.7) 
and the spatial representation found by the SWSH transforma¬ 
tion. Blue (dashed), L = 17/2, e r (i/;_); Green (dashed), 

L = 17/2, e r (i/) + )- Red (dashed), ” 0 ”: L = 33/2, e r {'tp_)-, Cyan 
(dashed), ”x”: L = 33/2, e r (ip+)’, Black (dashed), L = 65/2, 
e r (ijj-)', Magenta (dashed), ” 0 ”: L = 65/2, £ r ('i/j + ); Blue (solid), 
129/2, e r (^_); Green (solid), 129/2, £ r (^) + ). The thick 
black line corresponds to an expected convergence of 4th order in 
time. (See text for discussion) 



14 


F. BEYER, B. DASZUTA, AND J. FRAUENDIENER 


examples of the maximum absolute value of the 3-divergence V 0 j Q (expected to 
be 0 by Eq. (A. 2)) and probability conservation |1 — Q\ (again expected to be 
0) at a fixed number of time-steps (temporal discretization). This allows for the 
examination of the effect of varying the band-limit L on the evolution of the smooth 
initial data of Eq. (5.15). In Fig. 5(a) we observe that the numerical error is greatest 




(a) Maximum magnitude of 3-divergence 


(b) Probability conservation 


Figure 5. (Colour online) Test of the continuity equation and 
probability conservation for the case T = 1. We fix the time-step 5t 
to be 5/6400. N st indicates the number of time-steps taken from 
t = 0 (to the maximum t = 5). (a) Maximum over all sampled 
nodes in the spatial representation of |V a j“| (based on Eq. (A.2)) 

(b) Conservation of probability Q with time (based on Eq. (5.9)) 

In both sub-figures: Blue ”o”: L = 17/2; Red L = 33/2; 

Green L = 65/2; Black L = 129/2. (See text for 

discussion) 

for L = 129/2 and minimizes for L = 17/2 - the reason for this can be explained as 
follows: From the evolution equation (Eq. (5.6)) we see that coupling does not occur 
between distinct Z, m modes, only between fixed Z, m of ip±jm- Hence we see that if 
for given Z, m it is the case that ib-,lm\ t=0 = 4’+,im | t=0 = 0 then to within numerical 
tolerance these coefficients should remain 0 over the course of the evolution. Thus 
for a finite L we effectively provide a choice of initial data which can be viewed as 
exact at the specified L. The error arising from the spectral components of ip± can 
hence be entirely attributed to that of the SWSH transformation itself, which scales 
with increasing L as shown in Fig. 2. This behaviour can also be observed in Fig. 4 
by comparing the numerical error associated with components ip± at a fixed k for 
differing L. In Fig. 5 we examine the conservation of probability with the quantity 
11 — <51 where Q is calculated according to Eq. (5.9) and as T = 1 we may avoid 
rescaling to ip± as is required in more general cases. As no SWSH transformations 
are required to compute this quantity once ib±,im(t) is known we find that upon 
changing L the previous associated numerical error is not accumulated and the 
results for |1 — Q\ at differing L coincide. 

For a time independent bF we select the (finite L) initial data ib±.im(t)\ t=0 - 

ib—i _i = 2 , ib—i i = ib- a _3 = ib-s. a=l, 

' 2 ’ 2 2’2 ' 2 ’ 2 ' 2 ’ 2 

'lb— 7 _ 5 = i, 'll) _ 7 1 = —2i, '0— 1 7=1, 

T 2 ’ 2 ~ 2 ’ 2 2’2 


(5.16) 
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^+§,-§= 3 ’ V'+|,i=^+S,-|= 1 » ^+1,1=5, V’+|,-|=3, 

^+1,-1= 2, t;. 1 = 1. ^+5,1= 2 - ^+5,1 = 1, 

and J 7 we choose as: 

A,0 = 4, J 7 !,-! = —.Fi,! = -, ^4,2 = ^4,-2 = Yg- (5.18) 

Together with the above specifications we make use of auxiliary fields as defined 
in Eq. (5.4) with 4'±(t) = 0 => 4/± > i m (t) = 0 and the EOM decomposition of 
Eq. (5.5). In this case instead of comparing with an analytical solution we perform 
a self-consistent convergence based on Eq. (5.13). The result of this is presented in 
Fig. 6. We observe the expected 4th order convergence in time for the numerical 



k 

Figure 6. (Colour online) SCOT of numerical solution ip ± (spa¬ 
tial representation) for the case of static T (Eq. (5.18)). £ r (ip±) '■= 
e r (ip±(t; 5t k , i?, tp), ip±(t;5t k+1 ,’&,<p)) is displayed at t = 5 for a 
number of steps N = 2 k x 100. Blue (dashed), L = 11/2, 
e r ('i/'_); Green (dashed), L = 17/2, e r (ip+); Red (dashed), 

”o”: L = 33/2, £ r (^>_); Cyan (dashed), ”x”: L = 33/2, e r (ip + ); 

Black (dashed), L = 65/2, e r (ip_); Magenta (dashed), ”o”: 

L = 65/2, e r (ip- 1 -); Blue (solid), 129/2, e r (^»_); Green (solid), 

129/2, e r (ip + ). The thick black line corresponds to an ex¬ 
pected convergence of 4th order in time. (See text for discussion) 

solution at the band limits tested. In Fig. 7(a) and Fig. 7(b) we show 

representative examples of how the numerical solution obeys the continuity equation 
and probability conservation respectively 3 . Note that in order to calculate Q for this 
case we rescale fields as in Eq. (5.9) and use the full prescription of the PS method 
(see sec. 5.1). In contrast to the T = 1 case we now have a variable coefficient 
EOM and coupling between distinct l,m modes of i/j±jm occurs, thus sampling at 

3 The amount of raw data that can be generated at higher L can grow dramatically - during 

the calculation we thus only retain the coefficients V’i lm{t) a>t evenly interspersed points on the 

temporal grid as presented in figures. 
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a sufficiently high band-limit L in addition to choosing a sufficiently fine temporal 
grid is required in order to resolve conservation properties accurately. Finally we 



(a) Maximum magnitude of 3-divergence (b) Probability conservation 

Figure 7. (Colour online) Test of the continuity equation and 
probability conservation for the case of static T (Eq. (5.18)). We 
fix the time-step St to be 5/6400. N st indicates the number of time- 
steps taken from t = 0 (to the maximum t = 5). (a) Maximum over 
all sampled nodes in the spatial representation of |V a j Q | (based 
on Eq. (A.2)) (b) Conservation of probability Q with time (based 
on Eq. (5.9)) In both sub-figures: Blue ”o”: L = 17/2; Red 
L = 33/2; Green L = 65/2; Black L = 129/2. (See 
text for discussion) 

consider a time-dependent T. The selection we make is linear interpolation in 
time between initial < 7 ( 1 ?, y>) and final h{'d,ip) static deformations of § 2 . That is, 
:= (1 — + ( t/tf)h{d,tp) and tf = 5. The (real) functions g 

and h we choose to have non-zero gi^ m and hi^ m coefficients: 


to 

o 

o 

II 

0° 

5 

9 i,-i — ~g 1,1 — 2 ’ 

5 

H 1°^ 
' \ i-H ^ 

II 

(N 

II 

<N 

05 

(5.19) 

^o,o — 8, 

^i,-i = —/&i,i = — 

h/\ _2 — /14 2 — — 

4,2 1Q 


Initial data V , ±,im(i) | t _ 0 

is selected as in Eq. 

(5.16) and Eq. (5.17). 

Together 


with the above specifications we again make use of auxiliary fields as defined in 
Eq. (5.4) (note however that ’F-i-(t) ^ 0) and the EOM decomposition of Eq. (5.5). 
We again perform a self-consistent convergence based on Eq. (5.13). The result 
of this is presented in Fig. 8 . Once again, excellent agreement with the expected 
4th order convergence is observed with respect to the temporal scheme. We are 
again in a situation with a variable coefficient EOM, however now it also becomes 
non-autonomous. The continuity equation, together with probability conservation 
must be obeyed in this case also, and indeed we find that provided a sufficiently 
high band-limit L is chosen we may resolve these stated properties in our numerical 
solution (See Fig. 9(a) and Fig. 9(b)). 

5.3. Dirac equation: collapsing background geometry. Having analysed the 
numerical properties of our implementation of a pseudo-spectral method for the 
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Figure 8. (Colour online) SCCT of numerical solution i/>± (spa¬ 
tial representation) for the case of time-dependent T (Eq. (5.19)). 
s r (ip±) := e r {ip±(t w , St k , •&, ip), St k+1 , i?, yj)) is displayed at 
t = 5 for a number of steps N = 2 k x 200. Blue (dashed), 
L = 17/2, e T .(i/>_); Green (dashed), L = 17/2, e r (ip- 1 _); 
Red (dashed), ”o”: L = 33/2, e r (ip_)-, Cyan (dashed), ”x”: 
L = 33/2, e r {ip + )\ Black (dashed), L = 65/2, £ r (^>_); Ma¬ 
genta (dashed), ”o”: L = 65/2, £ r (^ + ); Blue (solid), 129/2, 
£ r {ip-)\ Green (solid), 129/2, £ r (^) + ). The thick black line 

corresponds to an expected convergence of 4th order in time. (See 
text for discussion) 


Dirac equation we now turn to potential applications. We consider a (2 + 1) di¬ 
mensional analogue of the usual Friedmann-Robertson-Walker (FRW) space-time 
in comoving coordinates. The particular physical scenario we wish to model is an 
imploding universe with pressure P and density p equal to zero. We take the scale 
factor to be a(t) = 1 — t = .F(t) -1 . Thus, we see that at t = 0 the background 
spatial geometry coincides with that of § 2 upon which initial data for the Dirac 
equation must be provided. 

For many practical purposes it is sufficient to consider an initial Gaussian state 
when discussing the dynamics of the Dirac equation. For a (n— l)-sphere embedded 
in K" the von Mises-Fisher distribution serves as an analogue to the planar Gaussian 
distribution [18]. The probability density function is given by: 

K n/ 2 -i 

Wx; x °' K) = p x) 
where k > 0, ||xo|| = 1 and l m {z) is the modified Bessel function of the first 
kind. The parameter n may be thought of as analogous to the reciprocal of the 
variance of the Gaussian distribution and x 0 as the mean direction about which the 
points x cluster. We choose initial data with average momentum 0 where j°| ( = 

/ 3 (x; x 0 | e=7r/ , 3 ^ =7r , 3 , 96) and j 1 = j 2 = 0. For the mass parameter we take p, = 1.2 
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0 0.5 1 1.5 2 2.5 


N st [xlO 4 ] 

(a) Maximum magnitude of 3-divergence 



0 0.5 1 1.5 2 2.5 


N st [xlO 4 ] 

(b) Probability conservation 


Figure 9. (Colour online) Test of the continuity equation 
and probability conservation for the case of time-dependent T 
(Eq. (5.19)). We fix the time-step St to be 5/6400. N st indi¬ 
cates the number of time-steps taken from t = 0 (to the maximum 
t = 5). (a) Maximum over all sampled nodes in the spatial rep¬ 
resentation of |V a j“| (based on Eq. (A.2)) (b) Conservation of 
probability Q with time (based on Eq. (5.9)) In both sub-figures: 
Blue ”o”: L = 17/2; Red L = 33/2; Green L = 65/2; 
Black L = 129/2. (See text for discussion) 


and evolve Eq.(5.2) and Eq.(5.3) on the temporal range t £ [0, 0.99]. As before we 
use explicit RK4 as the temporal integrator and upon performing SCCT (Fig. 10) 
find 4th order convergence in time as expected. Once again the continuity equation 
must hold and probability must be conserved. These properties are checked for 
the present numerical calculation in Fig. 11(a) and Fig. 11(b). For each chosen 
band-limit we observe growth in numerical error of several orders of magnitude as 
t 0.99. This is reasonable when we take into account that the spatial geometry is 
shrinking to a point (a(f) —> 0 as t —> 0.99). We now fix the band-limit as L = 129/2 
together with temporal grid N st = 6.4 x 10 4 . For the aforementioned resolutions 
snapshots of the dynamics are displayed on an Aitoff-Hammer projection [30] in 
Fig. 12. In Fig. 12(a) we see that the probability density j° is initially Gaussian 
in character and subsequently due to the rotational symmetry of the initial state 
about the mean direction Xo we observe evolution towards a ring-like structure 
indicating dynamics governed by dispersion (Fig. 12(b)). In Fig. 12(c) we see 
evolution towards the antipode of Xo, however we find that this does not coincide 
with a reconstruction of the initial condition in the limit t —> 1 (Fig. 12(d)). Observe 
that the average amplitude of lim t _>, 0 j°(t) increases dramatically; this is to be 
expected due to conservation of probability Q and the form of the volume element 
for this particular geometry used in its calculation (Eq.(A.ll)). 

6. Conclusion 

We have presented a new spectral algorithm for half-integer spin-weighted func¬ 
tions on S 2 based on spin-weighted spherical harmonics (SWSH) by extending the 
method for the integer case presented in [3,15]. Our implementation of this spectral 
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Figure 10. (Colour online) SCOT of numerical solution ip± 
(spatial representation) for the case of .F(f) -1 = a(t) = 1 — i. 
£ r (ip±) := e r {ip±(t w , 5t k , $, +), ip±(t; 5t k+1 , •&, ip)) is displayed at 
t = 0.99 for a number of steps N = 2 k x 4000. Blue (dashed), 
L = 17/2, e T .(i/i_); Green (dashed), L = 17/2, e r {ip .|_); 
Red (dashed), ”o”: L = 33/2, £ r (ip_)-, Cyan (dashed), ”x”: 
L = 33/2, e r {ip + )\ Black (dashed), L = 65/2, £ r (^)_); Ma¬ 
genta (dashed), ”o”: L = 65/2, £ r (i/j + )-, Blue (solid), 129/2, 
£ r {ip-)\ Green (solid), 129/2, e r {'ip + ). The thick black line 

corresponds to an expected convergence of 4th order in time. (See 
text for discussion) 


algorithm shows excellent agreement with the theory we have presented. Indeed 
we find that the numerical error scaling comparable to the integer SWSH algo¬ 
rithm [15] and that exponential convergence properties are displayed as anticipated. 
Furthermore, the expected algorithmic complexity of 0(L 3 ) is retained. In addi¬ 
tion, we have outlined how one may construct the 2 + 1 Dirac equation on a curved 
space-time, adapting the result to the 3-formalism and for a geometry with spatial 
topology of § 2 . Viewing the Dirac equation as an IVP we demonstrated how a 
pseudo-spectral approach can be combined with the half-integer SWSH algorithm 
and analysed several different situations with distinct spatial geometries, examples 
including linear interpolation between distinct spatial configurations with time and 
the case of an imploding 2 + 1 FRW-like model. The numerical solutions that we 
constructed obey anticipated convergence rates for the temporal solver used (ex¬ 
plicit RK4) and physical invariants (current continuity and probability) are shown 
to be preserved during the time-evolution to an excellent degree. 
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(a) Maximum magnitude of 3-divergence 


(b) Probability conservation 


Figure 11. (Colour online) Test of the continuity equation and 
probability conservation for the case of .F(t) _1 = a(t) = 1 — t. We 
fix the time-step St to be 0.99/6.4 x 10 4 . N st indicates the number 
of time-steps taken from t = 0 (to the maximum t = 0.99). (a) 
Maximum over all sampled nodes in the spatial representation of 
|Vcj Q |. (b) Conservation of probability Q with time. In both 
subfigures: Blue ”o”: L = 17/2; Red L = 33/2; Green 
L = 65/2; Black L = 129/2. (See text for discussion) 


Appendix A. The Dirac equation on a 3-dimensional Lorentz manifold 

Let V be a 77 ,-dimensional real vector space equipped with a quadratic form 77 with 
signature n— 2 , i.e., 77 can be represented in diagonal form as 77 = diag(l, — 1 ,..., — 1 ). 
Let c(V, 77 ) be the Clifford algebra associated to (¥, 77 ). This is the unique algebra 
with unit 1 and so called structure map 7 : V —>■ c(V, 77 ) such that for every v £ V 
the relation 7 ( 77 ) 7 ( 77 ) = r](v, 77)1 holds. Let (ei,..., e n ) be a basis of V and define 
7 i a b = rj(e ai eb) and j a = 7 (e a ). Then, by polarization of the defining relation, one 
finds the well-known Clifford-Dirac anti-commutator relations 

(7a, 76} = la'lb + Ibla = 2 Tj ab l. 

Irreducible representations S of the Clifford algebras have dimensions N = 2 "/ 2 
for n even and N = 2( n_1 )/ 2 for n odd. It is well known (see [26]) that these 
representations can be constructed recursively from the lower dimensional ones. 
The representation space S of the Clifford algebra c(V, 77 ) is called the spin space. 
Depending on the particular case, the spin space is equipped with certain invariant 
structures, see [26], [27], [31] for detailed discussions. 

The commutators of the generators 7 a b := 7 [ Q 7 &] = \ [la , 7 b] generate a Lie 
algebra with commutation relations 

[lablcd\ — Vcblad, Icalbd Vdblac ¥ I'Jdalbc- 

These are the commutation relations for the Lie algebra 0 ( 77 ) of the (pseudo-) or¬ 
thogonal group 0 ( 77 ) associated with the bilinear form 77 . This Lie algebra acts 
on the subspace 7 [V] C c(V, 77 ) as follows: for any skew bi-vector w ab and every 
vector v a we define 7 ( 777 ) := w ab "/ a b and 7 ( 77 ) = v a "/ a , then the action of 0 ( 77 ) is via 







A SPECTRAL METHOD FOR HALF-INTEGER SPIN FIELDS 


21 





Figure 12. (Colour online) Aitoff-Hammer projection of the 
probability density j° associated with the (2+1) dimensional Dirac 
equation on FRW background. Note that during the course of the 
time evolution the spatial radius a(t) of the spherical projection 
in each sub-figure is distinct, furthermore the colour scaling also 
varies. See text for discussion. 

commutation in c(V, rj) 

0 ( 77 ) x y[V] -A- 7 [V], (l(w),'Y{v)) fa 7 ( 111 ) 7 ( 1 ') — 7 ( 11 ) 7 ( 11 ;). 

We can see the corresponding action on V from the explicit calculation 
7 ( 111 ) 7 ( 1 ;) - 7 ( 11 ) 7 ( 10 ) = W ab V C {"fable - lelab ) = W ab V C ( l^bcla ~ Vaclb ) = 2 w a b V b ' y a . 

Thus, the action on 7 [V] corresponds to the linear mapping defined by 2w a b on V, 
which is clearly anti-symmetric with respect to 77 . The same Lie algebra also acts 
on the spin-space S via the representation of the Clifford algebra 

0 ( 17 ) x 5 A 5, ( 7 ( 10 ),^) fa 7 ( 10 ) 1 /;. 

In order to make use of spinors on the manifold .Jt we associate at every point x £ 
the Clifford algebra c (T x ^,g x ) and its representation space which we denote 
by S x . The rigorous construction is invariantly described in terms of (associated) 
vector bundles but we will not go into the details here. Instead we refer to [17] 
for a complete description 4 . The collection of all spin spaces S x forms the spin 
bundle S(^) over and sections of this bundle are referred to as spinor fields 
or simply spinors. Similarly, we obtain the Clifford bundle as the collection of all 
Clifford algebras c{T x ^,g x ). The structure maps defined at every point x yield 
a bundle map 7 from the tangent bundle Tb# to the Clifford bundle, assigning to 
every tangent vector v £ T x ^( an element "f x (v) in the Clifford algebra at x. 

The connection V on can be lifted to a connection also denoted by V on S{^//C) 
(see [37]). It can be uniquely characterized by the fact that the structure map 7 
and the invariant structures on the spin spaces are covariantly constant. With this 

4 Note, that in the general case there are global topological obstructions to the existence of 
these bundles, see [17,26]. However, since we are interested mainly in the 3-dimensional case 
where these obstructions do not exist we will not discuss them here. 
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connection one can define the Dirac operator on as follows: we pick a basis 
(ei, ..., e n ) of the tangent space T x ^ at each x £ . Now the anti-commutator 

relations are 

{7a, lb} = 2g ab l 

which hold at every x £ ^ with the appropriate definition of g a b and ~/ a . The 
Dirac operator S acting on spinor fields ip is defined by 

&!> ■= 5 ab 7(ea)V e > = 7 a V e >. (A.l) 


In this paper we assume the vector space V to have dimension 3 and to be 
equipped with a metric rj with signature The corresponding Clifford 

algebra c(V, 77 ) has a 4-dinrensional real representation [27,31], which may conve¬ 
niently be described within M( 2, C) as the real algebra generated by the ‘Dirac 
matrices’ 



'-1 O' 


' 0 1 


'0 i' 

7o = 

0 1 

, 7i = 

-1 0 

, 72 = 

i 0 


The spin space is C 2 (regarded as a 4-dimensional real vector space) and carries a 
sesquilinear form (Hermitean inner product) defined by 


(ipA) := ^*7o^, 

where ip* is the Hermitean conjugate of ip. For given spinors ip = 
this product is 

(lp,<P) = -Iplpl +7/>2^2- 

The generators are symmetric with respect to this product, i.e., we have 


and (j) = 


(tp, la4>) = ( 7 a,1pA}- 

Given two spinors <p and ip there is a naturally defined covector a = a(cp, ip) with 
components a a = (</>, 7aV’)■ When <p = ip, then this covector is 

a = [l^il 2 + I "021 2 , -1P11P2 - ^ 2 V’i,i(V’ 2 ^i - $1^2)] ■ 

Hence, it is real and null whenever (ip, ip) = 0. 

The commutators of the Dirac matrices q a are infinitesimal generators of the 
spin group Spin(l, 2). They can be written here in form 


7a7b - 7b7a = 2ie ab c 7 c . 


Choosing an orthonormal frame (eo,ei,e 2 ) on we can express the Dirac op¬ 
erator on with respect to the chosen frame explicitly. To this end we need the 
Ricci rotation coefficients for the frame defined by 


Y e a ^b — b a b e c . 

The condition that the structure map be covariantly constant translates into the 
equation 

r Q b C 7c = Tb7a ^ 7aTi, 

where the T^, which depend on the point x £ , are endomorphisms of the 

spin space S x . More precisely, they are linear combinations of the infinitesimal 
generators of Spin(l, 2). This equation can be solved uniquely for the T and yields 

Ta = -\T ba C e c bd ld - 

The connection defined in this way also leaves the complex structure and the 
sesquilinear product invariant. With these ‘spin coefficient’ matrices we can ex¬ 
press the Dirac operator acting on an arbitrary spinor ip as 

Slip = r] ab j a ( e b (ip) + T b ip). 
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The Dirac equation on can be obtained from a variational principle. For any 
given spinor field ip with compact support on SK we can write down the action 
functional 

2l[V>] := f i{ip,Sip) + niiii) d^f, 

where here and later on cL# denotes the invariant volume form on a manifold 
This defines a real number since the imaginary part can be rewritten as 

i f {ip, Si) + {ip, Shp) cL# = i / V a a a cL# = 0. 

J M 

Variation with respect to ip yields the Dirac equation 

i Sip + nip = 0. 

To every solution ip of the Dirac equation the covector a a defines a conserved 
current j a denoted by 

f = (Vb7 “V’)- 

Taking the divergence of this equation we find 

V a i a = (V a V>,7» + (V’,7 a V„V’) 

= {Sip, ip) + {ip, Sip) = {ini, Ip) + (i, ini’) = 0. 


(A.2) 


With this current we can define a conserved quantity in the usual way. Let t be 
a time coordinate for .S/, defined on an interval / = [0, T\ for some arbitrary T 
and let E be a 2-dinrensional manifold. We assume that I x E is embedded into 
as a 3-dimensional submanifold V via the embedding i : I x E 4 ^. Then 
it : E —> ,S/,, x i X ( t, x) embeds E as a space-like hyper-surface E t into M. Let 
S = <9E be the boundary of E, which i t embeds as a 1-dinrensional submanifold St 
into Now the boundary of V is V = E 0 U T U Et where T = U«e/ <5*. 

Integrating the divergence eciuation over V and using the generalised Stokes’ 
theorem we obtain 

o =[Va3 a dV=[ j a dV a , 

JV JdV 

where we denote the hyper-surface forms on the boundary dV by dV a - Introducing 
the future-pointing tinre-like normal t a to the hyper-surfaces E t and the outward 
normal n a to T we can write 


/ j a t a dT, = f j a t a dH + f j a n a dT. 

J So J S t J T” 

Defining the scalar “charge” (also interpreted as probability) Q and its flux J by 

Q{t)= f fta dE, J(t)= [ j a n a dS 

JY., JSt 


we can write the balance law 


Q{T) = 0(0) 


J (t) dt 


(A.3) 


(A.4) 


which holds for arbitrary values of T. With appropriate boundary conditions (for 
instance, when E has no boundary as we assume below) for the Dirac field one can 
make the flux vanish so that the balance law turns into the conservation equation 
for 0 

Q(T)=Q(0). 

We now specialize to the case ~ R x S 2 with the metric 

g = dt (Si dt — T~ 2 {t, d, <p) (dd S di? + sin 2 d dip S d^) , (A.5) 
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where ($, ip) are standard polar coordinates for the 2-sphere. The function T(t, 7 ?, ip) 
is a conformal factor relating the induced metric at every instant of time t on 8 2 
to the standard metric of the unit 2-sphere. We choose the orthonormal frame as 


eo = d t , e\=Td$ 1 e 2 = T cscfl d v . 


(A. 6 ) 


The non-vanishing Ricci rotation coefficients are 

ru° = r 22 ° = Tia^-i Tm 1 = - - 7- cottf). (A.7) 

J 7 J 7 sm i) 

From these we obtain the non-vanishing ‘spin coefficient’ matrices 


Ti 


1 

2 


Tip / sin •d 

. i Tt/T 


i T t /T 

—Tp! sini? ’ 


To = - 


T 7 # — Toot'd — T t /T 

J~t / T —J 7 # + J 7 cot i? 


(A; 8 ) 

In order to make use of the spin-weighted formalism we need to identify the spin- 
weights of the quantities involved. To this end we consider the infinitesimal rota¬ 
tions of tangent vectors to the sphere as induced by the (1, 2)-Lorentz group O(rj) 
at every point on jjt. Consider the infinitesimal rotation given by 2 = ^70 
acting as follows on the two frame vectors on the sphere 


7(ei) = 7i ^ 7) [7i2,7i] = 72, 7^2) = 72 H- - [712,72] = -71- 


In terms of the complex null-vector m 


= 7f ( ei “ ie2 ) is 


7 ( 771 ) 1 y i 7 ( 777 ), 


so ],7 i 2 generates the frame rotations m i-A e 1Q m. But it also acts on the spin-space 
at each point: for every spinor ip = [ipi, t/> 2 ] t we have 


1 

2 


712 ^ 


4>\ 

i 

-ip 1 

i>2_ 

= 2 

i> 2 _ 


Therefore, assigning spin-weight +1 to m results in the spin-weights —,) and +,) 
for ip\ and ip 2 7 respectively. Henceforth, we denote the spinor ip therefore by 7 p = 
ip+] T ■ Since J 7 is a scalar function it has spin-weight 0. 

Now we are in a position to write down the Dirac equation. Using (A.l), (A.6), 
(A.7) and (A.8) we find the two equations 


d t ip- = -ifiip- - \p'Tip. v - Td'ip + + ^ip~, (A.9) 

d t ip+ = ipip+ - ^dTip- - -T 7 + ^'0+- (A.10) 

The probability or charge density is 


j° = h4| 2 + KM 2 - 


We take V in the form (0 ,t) x § 2 with hyper-surfaces of constant time S T ~ § 2 
then the ‘charge’ integral as defined in (A. 3) yields the expression 


Qir) 



\ip + (r,d,p)\ 2 + |V>_(t,7?,V3)| 2 


sin 1 ) d$ dip 
T 2 (T,i9,(p)' 


Since E r ~ § 2 has no boundary, (A. 4) implies that for every r € (0, T) 


Q(t) = Q(0) = const. 


(ATI) 
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